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Stacks of intrinsic Josephson junctions in the resistive state can by efficiently synchronized by 
the internal cavity mode resonantly excited by the Josephson oscillations. We study the stability of 
dynamic coherent states near the resonance with respect to small perturbations. Three states are 
considered: the homogeneous and alternating-kink states in zero magnetic field and the homogeneous 
state in the magnetic field near the value corresponding to half flux quantum per junction. We found 
two possible instabilities related to the short-scale and long-scale perturbations. The homogeneous 
state in modulated junction is typically unstable with respect to the short-scale alternating phase 
deformations unless the Josephson current is completely suppressed in one half of the stack. The 
kink state is stable with respect to such deformations and homogeneous state in the magnetic field is 
only stable within a certain range of frequencies and fields. Stability with respect to the long-range 
deformations is controlled by resonance excitations of fast modes at finite wave vectors and typically 
leads to unstable range of the wave-vectors. This range shrinks with approaching the resonance and 
increasing the in-plane dissipation. As a consequence, in finite-height stacks the stability frequency 
range near the resonance increases with decreasing the height. 

PACS numbers: 74.50.+r, 85.25.Cp, 74.81.Fa 



I. INTRODUCTION 

Superconducting tunneling junctions are natural 
voltage-tunable sources of electromagnetic radiation due 
to the ac Josephson effect.^ As radiation from a single 
junction is very small, large-size arrays of artificially fab- 
ricated junctions have been used to enha nce power of 
electromagnetic radiation, see early reviews^' and more 
recent papersP^ The main challenge is to synchronize 
all junctions in the array. In this case the total emitted 
power is expected to be proportional to the square of the 
total number of junctions. 

Intrinsic Josephson junctions in the high- 
temperature layered superconducting materials^, 
such as Bi 2 Sr 2 CaCu208 (BSCCO), provide a very 
promising base for developing coherent generators of 
electromagnetic radiation which may operate in the 
terahertz frequency range. These materials have several 
important advantages in comparison with artificial 
structures made out of conventional superconductors 
including (i) the large packing density of the junctions, 
(ii) a large value of the superconducting gap (up to 60 
mev) which allows to bring the Josephson frequency into 
the terahertz range, and (iii) possibility to make very 
large arrays of practically identical junctions. 

The stack of junctions can be a powerful, coherent, 
and efficient generator only if the oscillations of the su- 
perconducting phases in all junctions are synchronized. 
Due to weak intrinsic interaction between the junctions, 
this is a challenging task. One possible way to synchro- 
nization is to use interactions between the junctions via 
the generated external radiatiorP. In this case, for ef- 
ficient coupling to the radiation field, a junction stack 
(mesa) must have small lateral size (< 10 /im) and con- 
tain a very large number of junctions (> 10 4 ). Such a 



mesa would be a frequency-tunable source with the con- 
siderable power conversion efficiency. The obvious tech- 
nological challenge of this design is the requirement to 
fabricate structures with such large number of almost 
identical junctions. This design has not yet been imple- 
mented in practice. 

A very promising route to efficient synchronization is 
to excite an inte rnal cavity resonance in finite-size sam- 
ples (mesas) P^l Being excited, the resonance mode can 
entrain oscillations in a very large number of junctions. 
The frequency of this mode is set by the lateral size of 
the mesa and for the resonance frequency in the tera- 
hertz range the width has to be rather large (40-100 p). 
The experimental demonstration^ and confirmation^^ 
of this mechanism marks a major advance in the quest 
for Josephson terahertz generators. 

In general, the structure of the coherent state and the 
mechanism of pumping energy into the cavity mode are 
nontrivial issues. Homogeneous phase oscillations at zero 
magnetic field do not couple to the Fiske modes. Such 
coupling can be facilitated by introducing an external 
modulation of the Josephson critical current densityP In 
this case the amplitudes of the generated standing wave 
and of the produced radiation are proportional to the 
strength of the modulation. 

An interesting alternative possibility has been demon- 
strated recent lypSHEH It was found that near the reso- 
nance an inhomogeneous synchronized state is formed. 
In this state, the stack spontaneously splits into two sub- 
systems with different phase-oscillation patterns, corre- 
sponding to the alternating phase kinks and antikinks 
statically located near the center. This leads to a static 
phase shift between the oscillations in the two subsys- 
tems varying from to 2n in a narrow region near the 
stack center. In spite of this c-axis inhomogeneity, the 
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oscillating electric and magnetic fields are almost homo- 
geneous in all the junctions. The formation of this state 
promotes efficient pumping of the energy into the cavity 
resonance. 

Another potential candidate for the coherent state pro- 
ducing strong emission is a homogeneous state in the 
external magnetic field, also known as a rectangular 
Josephson- vortex lattice. In spite of strong experimental 
efforts, the existence of this state has not yet been clearly 
demonstrated, except for small-size stacks at very small 
velocity.^ It was argued that in the large-size crystals 
the rectangular lattice is almost always unstable.^ For 
a finite-size system the stability analysis of the homoge- 
neous state has been performed recently^ and stability 
regions have been found. 

As large-size stacks have a huge number of degrees 
of freedom, stability of the coherent states with respect 
to small perturbations is an important and nontrivial 
issue. Linear stability analysis amounts to calculating 
the full frequency spectrum for small perturbations with 
respect to steady-state coherent solutions and verifying 
that there are no exponentially growing perturbations. 
The stability analysis allows to evaluate the range of pa- 
rameters where the coherent states are possible. For lin- 
ear arrays of point junctions the stability analysis has 
been performed in Refs. [T5] and [TIB In this case stabil- 
ity is strongly influenced by the external load. A large 
array of the small-size intrinsic Josephson junctions can 
be stabilized by the radiation field which acts similar to 
the external capacitive load. 7 

In this paper we perform a systematic comparative 
analysis of the linear stability of different coherent states 
in the array of extended Josephson junctions near the res- 
onance. We revealed two types of instabilities. The short 
wave-length instability corresponding to the alternating 
phase deformations develops for states which have re- 
gions of negative local time-averaged Josephson coupling. 
This instability is sensitive to the nature of the dynamic 
state. The homogeneous state in the modulated junctions 
is typically prone to this kind of instability. The long 
wave-length instability with wave lengths larger than the 
London penetration depth A appears due to the paramet- 
ric resonance excitation of the fast modes at finite wave 
vectors. Analysis of this instability is essentially identical 
for all dynamic states. The instability criterion depends 
on several factors including behavior of the resonance fre- 
quency shift with increasing the in-plane and c-axis wave 
vectors and the relation between damping of the uniform 
mode and modes at finite wave vectors. The most essen- 
tial parameters influencing stability include the shift of 
the Josephson frequency with respect to the resonance, 
the stack height, and the in-plane quasiparticle dissipa- 
tion. The paper is organized as follows. In Section |LT| we 
present the dynamic phase equations and coherent solu- 
tions near the internal resonance. In Sect ion [ill] we derive 
the linear dynamic equations for small perturbations with 
respect to steady states. In Section |IV| we consider the 
short-wave length instabilities of different steady states 



and present the numerical test for these instabilities in 
the stacks with modulated critical current. In Section 
|V]we describe the stability analysis with respect to the 
long- wave deformations and numerical verification of this 
analysis. The description of numerical simulations used 
to check some of our analytical results is presented in 
Appendix [5] 

II. PHASE DYNAMIC EQUATIONS AND 
COHERENT SOLUTIONS 

The dynamic equations for the Josephson-junction 
stacks have been derived in several papers^ and have 
been used in different forms in numerous simulation and 
theoretical studies PH These equations can be written in 
the reduced form as coupled equations for the phases ip n 
and dimensionless magnetic fields h„ = (h x>n ,h y>n ), 

d 2 if n , 2 , / d(p n - f \ ' fl , \ n 

-^-^-+(l-a\/ n ) \ v c-q — \-g\x) sixnp n — e ijz d t h^ n I =0, 

(1) 

£ 2 Vnhj, n - [l + Vab-^j [hj, n - e ijz diLp n ] =0. (2) 

Here i,j = x,y, eij Z is the Levi-Civita symbol (e xyz = 
—Syxz = !)• The units in these equations are selected as 
follows: 1/uip is the unit of time with io v being the plasma 
frequency, the c-axis London penetration depth A c is the 
unit of length, 4>o/(27rA c s) is the unit of the magnetic 
field. The function g(x) describes possible modulation of 
the Josephson current which was suggested as the way to 
couple to the internal resonanceP We consider both mod- 
ulated and unmodulated (g(x) = l) cases. The equations 
depend on four parameters, the layer-charging param- 
eter a, 22 two damping parameters, v c — 4ira c / (e c lu p ), 
v a b = ^o- a b/(s c j 2 ujp), and the ratio I = A/s, where a c 
and a a b are the components of quasiparticle conductiv- 
ity, A is the in-plane London penetration depth, and 7 
is the anisotropy factor. We will study a finite-size stack 
(mesa) containing N junctions with lateral sizes L x and 
L y . We consider the case L y L x . 

We tested some of our analytical results with numer- 
ical simulations. The numerical procedure is the same 
as in Ref. Q31 For completeness, details of the numerical 
simulations are described in Appendix |A"| 

A. Coherent states in zero magnetic field 

We consider the stack in the coherent resistive state, in 
which all junctions have identical voltage drops. In this 
state the dynamical phase can be written as^l 

ip n {x, t) « lot + 9 n (x,t). 

Further, we assume that the Josephson frequency ui is 
close to the in-phase resonance frequency Lo m = rmv/L x 



3 



and the resonance cavity mode is excited inside the mesa 
meaning that the oscillating phase has the large reso- 
nance contribution 6 n (x, r) ~ cos(mirx/L x ) cos(o;r + a). 
We will mostly focus on the experimentally relevant case 
of the fundamental mode m = 1. One can distinguish 
two particular cases: for a homogeneous (uniform) solu- 
tion the phases 6 n (x, r) are identical in all junctions and 
for an inhomogeneous (nonuniform) solution the phases 
9 n (x,r) vary from junction to junction. For the homo- 
geneous solutiorP, 9 n (x,r) = 8(x,t), Eqs. |l] ) and Q 
reduce to the sine-Gordon equation, 

In this case coupling to the resonance mode is only in- 
duced by the external modulation g(x). Representing 
the oscillating phase as 0(x,t) — Re[9 LU (x) exp(— iuir)], 
we obtain an equation for the complex amplitude 9 u (x), 



(lu + iv c uj) + 



if 



dx 2 



ig(x). 



(4) 



This equation has to be supplemented by the bound- 
ary conditions accounting for radiation. For symmetric 
mesas general boundary conditions can be presented in 
the following form 



dx 
dJL 
dx 



(L x ) = iC6 u (L x ) + i(9 u (0), 

(o) = -*C0 u (o)-*C0«(£x) 



(5a) 
(5b) 



The coefficients £ and £ depend on the particular geom- 
etry. For example, for the case of an isolated mesa on a 
metallic plate with thin metallic contact on the top 



u 2 L 



2i , 
.1 In 

2e c \ 7r 



u 2 L z 
' 2e c 



C 



[Jo(ku>L x ) + iN (k UJ L x )} , 



where k u = u/c, L z = Ns is the stack height, Jq(z) and 
No(z) are the Bessel functions, and C ~ 1. As £ and £ 
are small, near the resonance we can look for solution in 
the form 

8^ = ip cos(mirx / L x ) + v(x) 

where v(x) is the small correction accounting for the ra- 
diation boundary conditions. With this ansatz, Eq. Q 
becomes 

(uj 2 — uj 2 n +ii' c uj') ipcos (rmrx/L x ) 

12- \ 9 2 v . , 
+ + iv c uj) v + =ig{x) 

As v is small, we can approximately replace 



ui 2 v and look for correc- 



tion satisfying the mode-matching condition 



oj^v + d 2 v/dx 2 oc cos(mirx/L x ). In this case, v(x) can 
be found as 

v(x) « atp (x — L x /2) sm(rmrx / L x ) 

giving uj^v + d 2 v/dx 2 = 2aip cos{rrmx/L x ) (mn/L x ) . 
Substituting this result into the previous equation and 
taking projection to the mode, we obtain 



to 2 — ui 2 n + 2a (rmr/L x ) + iv c uj 



where g m = (2/L x ) J Q x cos(rmrx / L x )g(x)dx is the cou- 
pling parameter. The complex constant a has to be found 
from the boundary conditions ^ where in the righthand 
side we can neglect v{x). In this case the approximate 
boundary conditions for the correction become 



dv 



dx 
dv 
dx 



(L x )*i (-irc+c 



(o)«-i C + (-i) m C 



if, 



and they give identical results for a, 

a^(2i/rmr) k+(-l) m C 



The amplitude of the cavity mode can finally be repre- 
sented as 



where 



Re 



(1 + a r )uj 2 — w r 2 n + i (v r + v c ) u) 



(6) 



C + (-i) m C 



2luL 7 



[1 - {-l) m J {KL x )\ 



£cL x 

determines radiation contribution to the damping^H and 
4 



Loo 2 
ne c L x 



Im 

hi 



C + (-i) m C 
c 



ku,L 2 



+ (-l) m -N (k UJ L x ) 



determines the resonance frequency shift due to radia- 
tion. This small frequency shift is frequently neglected. 
However, it will be essential for the long-range stability 
analysis. Therefore, the homogeneous oscillating phase 
near the resonance can be represented as 



9(x,- 



Im 



g m exp(-iujT) 
(l+a r )u 2 — u 2 n +iv(jj 



mirx \ 
' ' /., J" 



(7) 

where v = v c + v r is the total mode-damping parameter™ 
The homogeneous solution is not the only possible co- 
herent state. A spectacular example of an inhomoge- 
neous coherent solution is the alternating-kink state re- 
cently reported in Refs. T2] and |T3] For this state the 
phase distribution is given by, 

9 n (x,r) = (-l) n e k (x)+8(x,T), 
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where, for the fundamental mode, the static-kink phase 
9k(x) changes from to n near the center. As the region 
of this change is extremely narrow^, in the equation for 
the homogeneous oscillating phase 9(x,t) one can ap- 
proximate 9k(x) with the step function 9f.(x) — > n Q(x~ 
L x /2), where 0(x) = for x < and Q(x) = 1 for x > 0. 
Within this approximation, such a state becomes equiva- 
lent to the stack with modulated current density with the 
modulation function gk(x) = sgn(cc — L x /2). In this case 
the homogeneous part of the oscillating phase is again 
given by Eq. 0, where for m = 1, g m -> g k ,i ~ 4/7T. 
While for the homogeneous solution coupling to the res- 
onance mode only exists due to the external modula- 
tion, for the alternating-kink state such coupling is self- 
generated. 



B. Homogeneous solution in magnetic field 



III. EQUATIONS FOR SMALL 
PERTURBATIONS WITH RESPECT TO 
COHERENT DYNAMIC STATE 

We will study the linear stability of the homogeneous 
solution with respect to small perturbations. A similar 
analysis has been done in a recent papei^ for the homo- 
geneous solution in magnetic field. We found, however, 
several important instabilities which were missed in this 
paper. 

We consider first the case of the homogeneous state 
in modulated junctions. With very good accuracy this 
analysis can also be applied to the kink state, because 
c-axis inhomogeneities in this state are located in a very 
narrow region near the center. Perturbing the homo- 
geneous solution, ip n (x, t) = lot + 9(x, r) + i? n (r, r) and 
h„(r, t) = h^ ' {x, r)+h„(r, r), we obtain the linear equa- 
tions for small perturbations #„(r, t) and h„(r,r), 



We also consider the homogeneous resonance solution 
induced by the external magnetic field h e applied along 
y-direction, 

(p(x, t) = LOT + h e x + 9(x, t) 

in a stack without the external modulation, g{x) = 1. 
The oscillating phase 9(x,t) — Re[0 w (x) exp(— icjt)] is 
determined by the following equation 



(c 2 + 



IV C UJ) 



d 2 9 u 
dx 2 



= i exp(— ih e x). 



We will mainly focus on the most interesting case of the 
fundamental mode and magnetic fields corresponding to 
a magnetic flux through each junction close to half flux 
quantum, h e — tt/L Xi providing the most efficient cou- 
pling to this resonance. The homogeneous solution for 
the junction stack is essentially identical to the corre- 
sponding solution for a single junctional, see also Refs. 
ITT1 and In particular, when the frequency is close 
to the resonance frequency o> m , the dominating contri- 
bution to the oscillating phase again has the resonance 
form given by Eq. ^ with 



5m ' gh.m 



L 



dx cos(mirx / L x ) e~xp(—ih e x) 



x Jo 



2i (1 - (-1)"' exp [-ih e L x ]) h e L x 
(h e L) 2 - (m7r) 2 



(8) 



being the coupling parameter due to the magnetic field. 
In contrast to the zero-field case, this parameter is a com- 
plex number. For the fundamental mode 



9h,i 



Mh e L x exp [-ih e L x /2] cos [h e L x /2] 
7T 2 - hlLl 



(9) 



In particular, g^ i — 1 for h e L x = n. This function is 
shown in the upper left plot of Fig. [I] 



d-d n 

I v r -Q--+g(x)C(x, T)0 n -eij Z dih jin 



0, (10) 



'dr 



hj, n ~ eijzdi'&n =0, (11) 



where C(x, r) = cos [wt + 9(x, r)]. Due to this oscillating 
cosine, perturbations are not monochromatic, oscillations 
with the frequency Q are coupled with the frequencies 
fi ± u>. At w > 1 we can look for solution in the form 

#„(r, t)r3 cos[q(n + 1/2)] cos (k y y) exp (— iflr) 

q,k y 

x ^ tf k ^(x)exp(ifiLJT), (12) 

/9=0,±1 

with q — Trm z /(N + 1), k y = irm y /L y 

and neglect other frequency components. For simple 
treatment of the c-axis modes, we assumed that the 
stack is bounded by metallic contacts which can be ap- 
proximated by ideal conductors.^ In this presentation 
fl = Q(k y ,q) is the complex eigenfrcquency which has to 
be found from equations ( 10 ) and (111 and the boundary 



conditions. The state is stable only if lm[Cl(k y , q)] < 
for all k y and for all g's from to tt. Separating the slow 
and fast parts in equations (10) and ( |11[ ), using relation 



g(x)C(x,T)-d n (r,T) 



2 ^ 

q,k y 



cos[q(n + 1/2)] cos (k y y) 



$k,o exp(j/3wr — £fir) + ^k.ff exp(— iVtr) 

/3=±1 p=±l 



and, excluding the magnetic fields, 

7 1 



l+£ 2 q 2 /[l-iv ab (Sl-Pu)] 
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with j3 = 0,±1 and q = 2sin(g/2), we obtain coupled 
equations for the slow and fast phase oscillations, see 
also Refs. [3 [HI and El 



nl g(x)C - tf M -r - 9 ,o^ 



G~ 



9{x) 



/9=±1 



0| 



r>- 2 ^ J5 , r -2 9 2 ^.js _ g{x)^ kfi 



(13) 

(14) 



where we introduced notations 



n 2 p = (n - ,3w) 2 /(l + cwf ) + iv c (O - /3w) , 

G^=i + ^ 2 /[i-^-^K(,]. 

Using 95 = ojt + Re exp(— zwr)], the time averaged 

cosine C{x) = (cosf) T can be evaluated as 



C{x) 



Im [0 u (a:)]/2 

ffm [(1 + a r )u 2 - cjgj cos(mirx/L x ) 
' [{l+a r )u 2 -ujU 2 + v 2 ^ 2 2 ' 



(15) 
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FIG. 1. (color online) Upper left plot shows the field de- 
pendence of the coupling parameter |<7h.i|, Eq. Right 
plot illustrates shapes of the local averaged Josephson cou- 
pling U(x)/U u , Eq. <|21[), for different magnetic fields and for 



Eqs. (pi) and (fTib have to be supplemented with the ^/K - w ) = 0.3. Lower /e/t plot shows blowup of the 



boundary conditions. At finite q, coupling to the external 
fields is negligible and with high accuracy we can use 
simple nonradiative boundary conditions 



region near the center for h e L x = ty to illustrate existence of 
the region with U(x) < which may lead to the short-scale 
instability. 



do 



dx 



for x = 0, L x 



(16) replacements 



The stability analysis reduces to computing spectrum of 
complex eigenfrequencies f2(k) from equations Eqs. (13) 
and ( 14 ) with nonradiative boundary conditions and find- 



ing out if there are regions in fc-space where lm[0(k)] > 0. 



9(x) 



(j=±l 0=±1 



exp [ij3h e x\ 



A. Case of homogenous state in external magnetic 

field 



in the righthand sides of these equations. The average in 
time cosine C(x) =(cosy) T can be evaluated as 



The above derivation can be directly extended to the 
case of the homogeneous state in external magnetic field. 
Perturbing the homogeneous solution, t p n (x , r) = lot + 
h e x + 9{x,r) + 'd n {Y,T), we obtain Eqs. ( fl0| ) and ( flT| ) for 
small perturbation $„(x, t), where now we have g(x) = 1 
and C(x, t) = cos (cut + h e x + 9(x, r)). Separating again 
the slow and fast components in the oscillating phase 
( 12 ), and using 



C(x, T)$ n (x, t) w i cos[g(n+l/2)] exp (— iOr) 
q 

x ^k,o ex P ( wr + h e x)] +^k,,9 exp (— ij3h e x) 
p=±i 

we obtain coupled equations for the slow and fast phase 



C(X): 



Re 



g h . m exp [ih e x] 



(1 + a r )u> 2 — uj 2 n + ivoj 



cos(mirx/L x ) 



oscillations, which are identical to Eqs. ( 13 ) and ( 14 1 with 



In the following sections we will analyze different insta- 
bilities of the dynamic coherent states. 



IV. SHORT WAVE-LENGTH INSTABILITY 

A. Modulated junction 

Consider first the region £q S> 1. In this case instabil- 
ity develops only for the homogeneous in the y-direction 
perturbations, k y = 0, and we only consider such per- 
turbations. In this regime the derivative term in the fast 
part ( 14 ) becomes small and can be neglected giving the 



G 



estimate 



'k,/3 



g(x) ( (n-fiuf 



1 



aq 2 



+ %v c (fi -PuS)\ d 



k,0, 



meaning that, roughly, i9k,/3 ~ ^k,o/w 2 <?C $k,o- Substi- 
tuting this estimate into the equation for the slow part 
(13), we conclude that coupling to the fast terms in this 
regime is weak and can be neglected. Therefore, the 
equation for the slow part becomes 



il 2 - U(x) 







k,0 



G 



,d 2 $ 



k.O 



9,0 Q x 2 



0, 



(17) 



where U(x) = g(x)C(x) determines the local plasma fre- 
quency. For the lowest mode we obtain 

u ( x ) = — 2 — ~ "^2 TV C0S (WAd)- ( 18 ) 

[w 2 - (1 + a r )u 2 ] +v 2 uj 2 

For large q we can use the simple boundary condition 
d$u.,o/d% — at the edges, x = 0, L x . As the lo- 
cal "spring constant" U(x) determines the local plasmon 
frequency fl\ oc oc ± y/U(x), the existence of regions with 
negative U(x) is a potential source of instability. This in- 
stability may be eliminated by the gradient term which 
we have to verify. Consider for definiteness the lowest 
mode and negative coupling parameter, g\ < 0. The 
simplest modulation giving such coupling is monotoni- 
cally increasing g{x). In this case at u> < u>\ the "spring 
constant" U(x) is negative at x < L x /2 and has min- 
imum near the edge x = where g(x) is the smallest. 
Therefore, the region for potential instability is located 
near the left edge, at i « in. As the most unstable 
mode is localized near x = 0, we can expand U (x) with 
respect to x 



g{x) cos{irx/L x ) « g(0) 



1 , ° 2 2 

1 + a\X —x 



with ai = </(0)M0), a 2 = (ir/L x ) 2 - g"(0)/g(0) 



giving the equation 



O 2 + U u ( 1 + aix - °^x 2 



k.O 



g: 



dx 2 



= 



with LL 



ff (0)| ffl |[c 2 -(l + a r V 2 ] /2 



[u 2 - {I + a r )uo 2 Y + u 2 u 2 
From this equation we estimate, 



> 0. 



n 2 



i 



2a 2 



C 



a 2 



UujG 2 



with 



«2 



_ (ir/Lf - g"(0)/ff(0) K~(l+«rV 2 ] +^ 2 
C/„G2" l + £ 2 q 2 /[l-tnu ab ] g(0)\ gi \[uj 2 -(l + a r )Lj 2 }/2 



and C ~ 1. Roughly, the system can only be stable 
if \a2U ul G~ 2 \ > 1 for all g's. As \G~ 2 \ has minimum 
at q = 7r = 2), it is sufficient to check the stability 
condition at this value of q. Skipping numerical factors 
on the order of unity, we can approximately write the 
stability criterion as 



9(0)91 (« 



(uj\ — UJ 2 ) + V 2 UJ 2 



< 



U 2 L 2 ' 



(19) 



This condition is very hard to satisfy because the right- 
hand side of this equation is very small due to £L X 3> 1. 
Therefore, the homogeneous solution is almost always un- 
stable with respect to alternating deformations localized 
near the edge with suppressed Josephson coupling. The 
formal reason for the large-g instability is that the lo- 
cal spring constant ( fl8| ) changes sign due to the factor 
cos(nx/ L x ). This instability would be completely elimi- 
nated if the modulation function g{x) would also changes 
sign in the middle. It is practically impossible to prepare 
such modulation artificially. The best artificial modu- 
lation for which the homogeneous state is stable with 
respect to the short-scale perturbations is steplike mod- 
ulation for which the Josephson current completely sup- 
pressed in one half of the stack, g{x) = for x < L x /2. 
In this case the plasma frequency is zero in this half and 



stability is achieved due to the gradient term in Eq. ( 17 ). 

We will demonstrate below that the alternating-kink 
solution is stable with respect to the large-g perturba- 
tion because it generates an effective modulation chang- 
ing sign in the middle of the stack making the spring 
constant positive in the whole stack. Another interesting 
case is the stack in the magnetic field corresponding to 
half flux quantum per junction. In this case the linearly 
growing contribution to the phase changes from to 7r 
across the junction which is similar to a sign-changing 
modulation. We will consider these cases in more detail 
in the following subsections. 



B. Absence of short wave-length instability for the 
alternating-kink state 

The stability analysis for the modulated system can be 
directly applied to alternating-kink state. It is crucial, 
however, that the effective modulation function gk{x) is 
close to a step function changing from 1 to —1 in the mid- 
dle of the stack, i.e., it changes sign in the middle which 
compensates change of sign in the factor cos(irx/L x ). 
This eliminates the short-scale instability. Indeed, the 
effective "spring constant" (18), which determines the 



local oscillation frequency, can be evaluated as 



U k (x) = g k (x)C(x) 

_ gfc.l K - (1 + Qr)^ 2 ] /2 

~ [lo 2 - {I + a r )u 2 ] 2 + u 2 u 2 



cos(ttx/L x )\. (20) 



It only touches zero in the midpoint and never goes nega- 
tive within the stack. Therefore, in contrast to the homo- 



7 



N=50,L =25^„«=100 ' 

1 X J' 


|F )"1 1 

1. 1. 


v c =0.01 ,v ab =0.2 

kink state 


/■ 


■ i r 


stable uniform ^ 




. (steplike modulation)— 


4 1 ■ 


unstable uniform A — 




"(linear modulation)^<~^ s 















y 







0.02 



12.28 12.36 12.44 12.52 

V 



50 100 150 200 



t = 


t = 20T, 


t = 40Tj 


































FIG. 2. (color online) Upper left plot shows current-voltage 
dependences near the resonance for three dynamic states: 
(i) unstable uniform state in the stack with linear modula- 
tion of the Josephson current, jj(x) = jjo(l + r(2x/L — 1)) 
with r = 0.4, (ii) kink state, and (iii) stable uniform state in 
the stack with steplike modulation of the Josephson current, 
jj(x) = jjoQ(2x/L — 1). Dashed lines show corresponding 
theoretical curves. Instability for the first state was triggered 
at the point j = 0.13 marked at the plot. Lower plots shows 
snapshots of the phase distributions in the eight bottom junc- 
tions for t = 0, 20Tj, and 40Tj, where Tj = 2-k/uj w 0.5 is 
the period of Josephson oscillations for the uniform state. All 
phases are shifted to the range (— 7r,7r). The snapshots il- 
lustrate development of instability near the left edge. Upper 
right plot shows the time evolutions of the difference between 
phases in the 6th and 5th junctions at the left edge. This 
difference changes from corresponding to the uniform state 
to 2-7T corresponding to the kink state. 



geneous state in modulated junction, the alternating-kink 
state is stable with respect to short-scale perturbations. 



FIG. 3. (color online) Stability range of the homogeneous 
state in magnetic field with respect to the short-scale phase 
deformations based on Eq. ( 22 1 with the listed representative 
parameters. 



ops near the edge with suppressed Josephson current, as 
the analytical analysis predicts. The upper right plot in 
Fig. [2] shows the time evolutions of the difference between 
phases in the 6th and 5th junctions at the left edge. We 
can see that this difference evolves from corresponding 
to the uniform state to the value 2n corresponding to 
the kink state. This transition is accompanied by large- 
amplitude oscillations with the period much longer than 
the period of Josephson oscillations. These oscillations 
are a consequence of the small c-axis dissipation param- 
eter, v c — 0.01, which we used in our calculations. 

We also verified that the uniform state remains stable 
in the stack with steplike modulation of the Josephson 
current, jj(x) = jjo for x > L x /2 and jj(x) — for 
x < L x /2. The current- voltage dependence for this state 
is also shown in Fig. [2] 



Test of short-scale instability in modulated 
stack with numerical simulations 



D. Homogeneous state in magnetic field 



We verified the short scale instability in modulated 
junctions using numerical simulations described in Ap- 
pendix [X] First, we probed the stability of the uniform 
state in the stack with linear modulation of the Joseph- 
son current, jj(x) = j,/o(l + r(2x/L — 1)). We used 
the modulation parameter r = 0.4. If we start from 
the uniform n-independent state and solve the dynamic 
equations without noise than instability does not develop 
and we can trace the current- volt age dependence corre- 
sponding to this uniform state, see Fig. [2j However, if 
we add to the phases a small alternating perturbations, 
(—l) n S(p, than we observe that the uniform state blows 
up and, after extended time evolution, it converges to 
the dynamic kink state. The initial stage of this time 
evolution is illustrated in the lower part of Fig. [2j From 
the phase snapshots we can see that the instability devel- 



Due to the complex coupling function in the case of fi- 
nite magnetic field, the short-range stability has features 
which are special for this case. In the regime £q 1 
the coupling to the fast phase can be neglected, as for 
the modulated stack case. This means that the slow 



part again obeys Eq. (171 in which the local "spring con- 
stant" U(x) now is simply given by the average cosine, 
U(x) = C(x). For the fundamental mode, we obtain 



U(x) = U u 



— sin 



with LL 



h 
2 

2h e L x cos [h e L x /2] 



2 



(21) 



{h e L x y ( w a - uir + 



For simplicity, we omit here the factor (1 + a r ) in front 
of u> 2 which has very little influence on the short wave- 
length stability. Shapes of this function at different 
fields is illustrated in the upper right plot of Fig. [T] for 
vuj/[uj\ — uj 2 ) — 0.3. Note that U(x) always changes sign 
at the center, x — L x /2, see, e.g., lower left plot in Fig. 
[I] However, for vuj <C uj\ — uj 2 the region of negative 
U(x) is very narrow and its existence does not automat- 
ically imply instability. We analyze the central region 
as the most prone to instability. For uj 2 — lu 2 ^> vuj 
the "spring constant" behaves near x = L x /2 as U(x) 



Eq. (fl7|) becomes 



h e x\ irx/L x with x = x — L x /2, and 



n 2 



VUJ \ 7TX 

2 2+h e x — 

i-uj 2 ) L x 



7. r - 2 ^k,0 _ 



From this linear-oscillator-type equation, we derive equa- 
tion for the complex cigenfrcquency 



n 2 /(l + aq 2 ) + iv c n 



v 2 uj 2 tt/L x 
ih e (uj 2 - uj 2 ' 



lirh e /L x 



which gives the stability criterion 



nh e /L x 



> 



V UJ 7T / L x 



Y 2* 2 (1- cos q)U u Ah e {uj 2 ~uj 2 f 

For the most "dangerous" mode at q = ir, assuming uj\ - 
uj 2 3> vuj, we obtain the following stability criterion 



uj 2 > 



n\g lhl \e 2 v W 

mi*, 



1/5 



(22) 



Large value of I 2 in the righthand side is compensated 
by small value of z/ 4 u; 4 . Note that increasing dissipation 
reduces the stability range. Taking typical values h e ss 5, 
uj\ rj 5, h e L = 7r, £ = 150, and v — 0.002, this inequality 
gives (uji — uj) /uji > 1.3 • 10~ 3 and for larger dissipation, 
v = 0.01, (uj\ — uj)/uj\ > 0.005. This conditions are not 
too restrictive. A representative stability region in the 
magnetic field is illustrated in Fig. [3| We conclude that 
the homogeneous state in the magnetic field correspond- 
ing to half flux quantum per junction remains stable with 
respect to the short-scale deformations if the frequency 
is not too close to the resonance. 



V. LONG- WAVE LENGTH STABILITY 

We analyze now the acoustic type instability at very 
small k y and q. With very minor modifications, this anal- 
ysis applies to all dynamic states considered in this pa- 
per. Obviously, it is most relevant for the systems which 
are stable with respect to the short-length deformations. 
Consider equation ( 14) for the fast components $k,±i- At 



small k y and q this equation gives the resonance solution. 
To obtain approximate solution for i?k,±ij we will keep 
only the resonance term, $k,/3(x) ~ V'k,/? cos(rmrx / L x ) . 
Further analysis shows that the coordinate dependence 
of $k,o is weak and can be neglected. In this case, the 
mode amplitude tpu.,/3 can be found following the same 
reasoning as for the homogeneous solution leading to 



V'k,/? 



gm^k,o/2 



o 2 



Using this result, we present Eq. ([13]) for the slow com- 
ponent $k,o as 



h 2 -g{x)C{x 



n~ 2 k 2 



GZ 



.d 2 d K0 



9{x) 



k,0^-,,0 g x 2 

g m $k,ocos(mirx/L x ) 



^ Q 2 

/3=±1 "(9 



■*2) 



(23) 



The typical length scale for its variation, Iq 

i 



o2 _ r- 2 h? 



G 



9,0 



exceeds the stack width L x 



because we consider the case of small k y , q, and f2 when 
Gqfi ~ 1 and il,k y <C l/L x . This allows us to ne- 
glect x dependence of #k,o- In this case g(x)C(x) and 
g(x) cos(mirx/L x ) can be replaced by their averages over 
x 



{9{x)C{x)) x 



dig [(1 + CXy)" 2 ~ /4 

[(l+a r )uj 2 -uj 2 n } 2 + (v c + v r ) 2 uj 2 



(g(x) cos(mnx/L x )} x = g m /2. 

We also can neglect the charging effects because at the 
typical wave vector q ~ tt/N and a < 0.1 the charging 
correction aq 2 is tiny and has only minor influence on 
stability criteria. In these approximations, the equation 
for Cl(q) becomes 



n 2 - 



i 



, 2 5 -o 1 + i" c n = fi^Q{^ k) 

k z z /[l-i\lv ab \ 8uj m 



(24) 



with k z = 2£sm(q/2) w £q, k = (k y , k z ) and 
2uj m [(1 + a r )uj 2 - uj 2 m ] 



G(n,k) = - 



[(1 + a r )uj 2 — uj 2 ^] + v 2 uj 2 



f)=±i 



(n-Pujy 



k 2 



k 2 z 



1- 



-+iv c (Cl—fioj) 



If we eliminate the radiaton corrections, v r and a r , which 
only appear at very small k y and k z then this function 
satisfies the translational invariance condition C? (0, 0) = 
0. Note that we can not take true limit k y ,k z — > 0, 
because we consider a finite-size stack with geometrical 
sizes smaller than the wave length of outside radiation. 

As the Josephson frequency uj is close to the cavity- 
mode frequency uj m , we can keep only the dominating 
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resonance terms. To simplify presentation of the second 
term, we introduce the shift of the plasma frequency at 
finite k y and k z with respect to the homogeneous mode, 



A k = uj p (k m , 0, 0) - u p (k 
parameter, 



k y ,k z ), and mode damping 



1 



A , ^ 



UJ m — 



ky 



1 +**/(! 



ky 



iv ab (0 - /3w)) 



1 + ">m 



2 1 



Introducing also the resonance detuning 6 U 
we simplify (SI, k) as 



A, 



(O + iv k /2Y 



(25) 



with <5 r = a r u m /2 being the radiation shift of the reso- 
nance frequency and 



i v abk\ 



1 



(26) 



is the damping of the plasma mode at finite wave vector 



(we neglect small terms on the order of k z k y ) 

Due to the resonance structure of 5(51, k), in the limit 
Ak, \S^\ <C gm/V^m h typically exceeds the righthand 



side of Eq. ( 24 ) and the dispersion equation is approxi- 
mately given by (S7, k) = 0. In this approximation we 
obtain the following result for the eigenfrequencies 



n±(k) 




(27) 



where we introduced new notations S u = S u + S r and 
Ak = Ak — S r to absorb the radiation frequency shift. 
Only the mode S2 + (k) is potentially unstable. The insta- 
bility takes place when the expression under the square 
root is negative and the imaginary square root exceeds 
the first term. This leads to the following condition for 
the instability in the (k y , k z ) plane 



(S„ + A k )(£-A 



> 



(28) 



Analysis of this criterion shows that the long-range sta- 
bility is determined by several factors including behav- 
ior of the plasmon frequency shift Ak and the relation 
between the damping of the homogeneous mode v and 
damping of excited modes at finite wave vectors, v k . As 
the shift Ak may take both positive and negative values 
depending on k y and k z , it is more transparent to find 
the instability range for this parameter, 



A, 




(29) 
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ol=0.0015, v,=0.02 
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FIG. 4. (color online) (a) Evolution of the instability region for 
the long-wave deformations in the k z -k y plane with approach- 
ing the resonance from below. Unstable regions are inside the 
domes. Used representative parameters are shown in the plot. 
The curves are marked by the value (uji — uj)/uj 1 = — 5 u /cj m . 
Dashed line shows dependence k y /ki = k z /^l + u^ b ujf cor- 
responding to condition to p (k 1 ,k y ,k z ) = cj p (fci,0,0). The 
points illustrate discrete wave vectors for a finite-size stack. In 
this example, the system becomes unstable for (u m — uo)/uJm 
between 0.04 and 0.08. (b)Evolution of the instability region 
in the k z -k y plane at fixed frequency with increasing the in- 
plane dissipation parameter v a \,. 



In particular, there is no instability in the fc z -region sat- 
isfying the condition 



Vk > \Su. 



4|<L 



which explicitly can be written as condition for k z 

k. 



1 + 2 



(30) 



4|S M | 



Qualitatively, we can conclude that the in-plane dissipa- 
tion suppresses instability until v a i,uj m < 1, while the 
radiation damping enhances the instability. As for a 
finite-size stack the fc z -values are limited from below, 
k z > £ir/N, the above equation also determines the sta- 
ble frequency range for a stack with given size. We can 
see that this range expands with decreasing N. 

In small-fc z range where the condition (30) is not sat- 
isfied there may be a range of k v w here the system is 



unstable. We can find from Eq. (29) an explicit presen- 
tation for this range 



A 
OJr, 



4:S U UJ r< 



< 



\ 




(31) 



The instability regions in the (k y , fc z )-plane for the dif- 
ferent Josephson frequencies near the resonance are il- 
lustrated in Fig. [4^,. One can see that the instability 
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FIG. 5. (color online) Stability range with respect to the long- 
wave deformations near the resonance as a function of the 
stack height based on Eq. ( |33[ ). Two regimes are illustrated 
in plots (a) and (b), which are controlled by the in- plane 
dissipation. Dashed lines show stability boundaries without 
radiation corrections. In the calculations we used the follow- 
ing representative parameters, v c — 0.005, oj\ = 10, t = 130. 
For radiation parameters we used a r = 3.3 ■ 10~ 4 L 2 /A a 5 and 
Ur/uJi = 5.2 ■ 10~ 4 L Z / \ a b, where the numerical coefficients 
were estimated assuming A a b = 0.25/im, L x = 80/im, e c = 12, 
and hx(C/k 



region rapidly shrinks with approaching the resonance. 
Vanishing of instability at large k z is caused by increas- 
ing mode damping due to the in-plane dissipation. This 
is illustrated in Fig. |4Jd where we plot the instability re- 
gion at fixed Josephson frequency for different in-plane 
dissipation parameter v a b- We can see that the insta- 
bility regions shrinks with increasing is a f,. Existence of 
the instability region in the (k y , fc z )-plane does not au- 
tomatically imply instability in real junction stacks. In 
a finite-size stack a discrete set of the wave vectors is 



allowed, k Vtm — mn/Ly, k z ^ n = nlir/N. The system is 
only unstable if at least one of discrete pairs (fc y ,m, k Z)n ) 
falls inside the instability region at given frequency, as 
illustrated in Fig. [4ji. 

At finite radiation corrections and at sufficiently strong 
in-plane dissipation the instability region may vanish 
completely in some frequency range. The condition for 
absence of instability can be written as 



max 

k z 



1 + "lb UJ ?n 



V 

4<L OJr, 



A 




< 



Finding the maximum leads to the following range of 
frequency at which there is no instability at all 




(32) 



with W a b — \A + < jJ m v tb JrijJ 'niVab- This region only exists 
if inequality Lu^ n v a i,a r > v r is satisfied, leading to the 
following condition for the in-plane dissipation Lo m v a b > 
2n/\n[C/k LJ L z }. 

The largest value of k z at which the instability bound- 
ary crosses k v = 0, k z g, can be found as 



"z0" 



-v c v ah 



1 + (u> m VabY 



hw m i/ ol a r \--~ 



(33) 



A finite-height stack is stable with respect to the long- 
range deformations if the minimum wave vector k z m i n = 
ir£/N = ir\ a b/L z exceeds the maximum between the two 
frequency-dependent wave vectors k Zyi and k z o defined by 



Eqs. (30) and (33) leading to the following criterion 



N < 



max(fc 2i „fc z0 ) 



(34) 



This condition gives to the frequency-height stability di- 
agrams illustrated in Fig. [5] As the radiation correc- 
tions a r and v r are proportional to the stack height, 
these boundaries have to be computed self-consistently. 
A simple estimate for the wave vector k z0 can be ob- 
tained for weak in-plane dissipation u) m v a b -C 1 away 
from the resonance at 5 U < 0, \5 U \ S> v, oj m a r . In this 
case k z0 « 2|<5 W | meaning that the finite-height becomes 



unstable at \5 u \/w m > (tt£/A) 2 /2. 



In presence of radiation corrections, two regimes exists 
depending on strength of the in-plane dissipation. At 
small uj m v a b there is always the critical stack height above 
which the system becomes unstable. Coupling to the ra- 
diation decreases this critical height, i.e., it enhances the 
instability. This is in contrast to the small-size stacks 
away from resonances^, where the coupling with outside 
radiation stabilizes the synchronized state. This behavior 
changes at large inplane dissipation. In this case a range 
of frequencies exists within which the system remains sta- 
ble for all stack heights. Also, in this regime coupling 
to radiation somewhat increases the critical stack height 
away from the resonance. 
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FIG. 6. Probing the long-range instability using numerical 
simulations of the coarse-grained model, (a) The current- 
voltage dependences for three sets of parameters below the 
resonance voltage V\ = 4ir. The horizontal arrows mark on- 
sets of instabilities in simulations and the vertical bars mark 
the theoretical estimates described in the text. The insta- 
bility leads to appearance of a small additional bump in the 
current- voltage dependence, (b) The time evolution of the dif- 
ference between the bottom and top phases at the left edge. 
One can see that for j > 0.121 this difference decays with 
time indicating stability of the homogeneous state, while for 
j < 0.120 this difference remains finite. The inset shows the 
blowup plot of the long-time decay of the phase differences 
for stable states. The stripelike appearance of the curves is 
due to the rapid oscillations with the Josephson frequency. 



A. Test of the long-range instabilities with 
numerical simulations 

We made several approximations in our analytical 
derivations, which allowed us to obtain relatively sim- 
ple criteria for the long-range stability. To verify va- 
lidity of these approximations, we checked some of the 
analytical results numerically. As suggested by the rep- 
resentative phase diagrams shown in Fig. [5j the long- 
range instabilities are only expected in rather tall stacks 
with heights N exceeding 1000 junctions. It is very diffi- 
cult to simulate such tall stacks directly. Fortunately, if 
we only interested in the long-range instabilities, this is 



not necessary. For this purpose, we simulated a coarse- 
grained model with the step in c-direction Sz contain- 
ing many junctions. As demonstrated in the Appendix 
|A"1 by change of variables, this coarse-grained model can 
be reduced to the original model with reduced parame- 
ter t, i — > £ cg — £/Sz. This trick allows us to use the 
same code with different parameters to probe the long- 
range stability of very tall stacks. In this case, the num- 
ber of layers in the model N is replaced by the number 
of numerical slices -/V cg = N/Sz = N£ cg /£. In numer- 
ics we use the two-dimensional model which only allows 
us to check our analytical results in the simplest situa- 
tions when the instability is homogeneous in y-direction 
(k y = 0). We also neglected the layer-charging effect, 
a = 0, and did not take into account the radiation cor- 
rections, a r — i> r = 0. Having in mind to probe the 
long-range stability of the kink state, we use the modu- 
lation function g(u) = sign(u — L x /2). With such mod- 
ulation function the system is stable with respect to the 
short-scale perturbations 

To probe the long-range stability, we numerically 
solved the dynamics equations for increasing transport 
current in the voltage range corresponding to the Joseph- 
son frequencies close and below the fundamental reso- 
nance. We added small deformation Sip(u, n) cx cos(7r(n— 
1/2) /N) at the beginning of every run for new value of 
the current and monitored the time evolution of the dif- 
ference <p(0, l)—tp(0, N cg ). Figure|6]shows results of simu- 
lations of the coarse-grain model, which reveal the long- 
range instability. Figure [6^ shows the current-voltage 
dependences near the resonance for three sets of param- 
eters: (i) £ cg — 5, L x — 1.25Xj8z, and N cg = 50, (ii) 
£ cg = 10, L x = 2.5XjSz, and A cg = 50, and (hi) £ cg = 5, 
L x = 1.2bXjSz, and A cg = 40. In all cases we used the 
dissipation parameters v c = 0.01 and u a b = 0.2. The 
values of the total stack heights N — N cg £/£ cg listed in 
the legend were obtained assuming £ = 150. For all three 
cases we observe the long-range instability which leads 
to appearance of a small additional bump in the current- 
voltage dependence. Development of the instability is 
illustrated in Fig. [6]d in which we show time evolution 
of (f(0, 1) — ip(0, N cg ) for different currents for the first 
set of parameters. We can see that for for j > 0.121 
the initial perturbation decays with time indicating sta- 
bility of the homogeneous state, while for j < 0.120 the 
perturbation does not decay. The instability onsets are 
marked by the horizontal arrows in Fig. [B£i. The sets of 
parameters (i) and (ii) correspond to simulations of the 
same physical system with two different coarse-graining 
parameters, Sz — 30 and 15 (assuming £ — 150). We ob- 
served that in both cases the instability develops at the 
same voltage indicating that the coarse-graining does not 
influence much the long-range stability. The set (iii) has 
the same parameters as the set (i), except for the smaller 
height. We see that the instability moved to the lower 
voltage, as expected. 

We now compare the location of the instability onset 
with the analytical predictions. At k y = the condition 
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FIG. 7. Snapshots of the distribution of the electric field for the steady state in the unstable region for the first set of parameters 
in the previous figure at j — 0.12 over the half period of the Josephson oscillations. One can see the contribution from the 
nonuniform mode. 



of stability is given by ir£/N > k z o and, without the ra- 
diation corrections, the formula (33 1 for k Z Q significantly 
simplifies, 



\M 



v c v ah I , for a r , v r = 0. (35) 



Moreover, for the parameters we used in simulations the 
terms with v c are negligible and with high accuracy we 
can estimate the shift from the resonance where the in- 
stability is expected as 

|<L| « ( Wl /2)(vrf/A0 2 - (wi/2)(^ cg /^Vcg) 2 . 

This gives |<5J w 0.62, V inst = ^ - \S U \ « 11.95 for the 
first and second sets of parameters and \6 U \~ 0.97, Vinst ~ 
11.6 for the third set. These values are shown by the 
vertical bars in Fig. [6^,. We see that in the simulations 
the instabilities appear exactly where they are predicted 
analytically. This gives us a confidence that the used 
approximations are legitimate. 

With simulations we can go beyond finding the loca- 
tion of the instability onset. We can also find the finite 
dynamic state after the instability develops. To under- 
stand the structure of this finite state, we present in Fig. 
[7] snapshots of the distribution of the electric field for 
the first set of parameters at j = 0.12. Analyzing these 
snapshots, we conclude that the instability leads to the 
state in which the oscillating phase is a superposition of 
two modes, 

<WO) = [#i,o + i#i,i cos(7r(n-l/2)/iV)] cos(ttx/L x ), 

where the amplitude of the nonuniform mode, #i,i, con- 
tinuously grows starting from zero at the instability on- 
set. 



VI. SUMMARY 

In conclusion, we found that dynamical states syn- 
chronized by the internal cavity resonance are prone to 
two very different instabilities. The short wave-length 
instability develops for states which have regions of neg- 
ative time-averaged Josephson coupling. In particular, 



the homogeneous state in stacks with modulated Joseph- 
son coupling typically has this type of instability. The 
homogeneous state in the external magnetic field H < 
<S>o/(sL x ) has this type of instability close to the reso- 
nance and the instability range widens with decreasing 
field. The long wave-length instability appears due to 
the parametric resonance excitation of the fast modes at 
finite wave vectors. The instability criterion depends on 
the relation between the damping of the homogeneous 
mode and modes at finite wave vectors. Finite-height 
stacks are stable sufficiently close to the resonance. The 
instability region typically shrinks with increasing in- 
plane dissipation. 
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Appendix A: Numerical simulations and 
coarse-graining procedure 

For numerical simulations it is convenient to present 
the dynamic equations in the form of the time-evolution 
equations for the reduced c-axis electric fields (e„), 
phases (^ n )> in-plane supermomenta (fc n ), and magnetic 
fields (h n ), 



v a b 



de n 
dr 

dr 
dk n 

~3t 



= -v c e n - g(u) sincpn + 



dhn 
du 



— [k n + h n h n —i] , 
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(A2) 
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The units in these equations are different from units 
used for analytical calculations: unit of length is the 
Josephson length Aj, unit of supermomentum is 1/A,/, 
unit of magnetic field is $0/(2^7 A 2 ), and unit of elec- 
tric field is <5>qlu p / (2-kcs). All parameters are assumed to 
be y- independent. Therefore we only probe instabilities 
uniform in this direction. We also neglected the layer- 
charging effect, a — 0. Above equations are solved for 
stack containing N junctions with < u < L x assuming 
simple non-radiative boundary conditions at the edges, 
k n = 0, d(p n /du — ^fI/2£ 2 at u = 0, L x , where I = jL x 
is the total transport current 



FIG. 8. Illustration of the staggeredgridused for numerical Solution of these equations is implemented using the 

solution of the dynamic equations jAl)-||A4}. following implicit numerical scheme: 



• Space and time discretizations are performed using a staggered grid. For coordinate, ip n (u,r) and e n (u,r) are 
defined at the points u — (j — l/2)d u , while k n {u,r) and h n (u,T) are defined at the points u = jd U: see Fig. [8] 
For time, e n is defined at r = md T while (p n , k n , and h n are defined at r = (to + l/2)d T . 

<r /2 = P« [C? " V2)du, (m + l/2)d T ] , e%j = e n [(j - l/2)d u ,md T ] 1< j < J 

Cf /2 = k n [(j - l)du, (m + l/2)rf T ] , Ct V2 = [(j - IK, (to + l/2)d T ] , K j < J + 1 



We discretize equations as 



m+l _ m 

d T 

m+3/2 m+1/2 



d T 



i+3/2 



-1/2 



n,3 



i+3/2 



m+l , m 
n.i ~ n,j 



m+l 



, m+1/2 , m+1/2 
m+1/2 n,j + l ~ 



n,3 



1 

Vab 



, m+3/2 , m+1/2 , m+3/2 , m+1/2 , m+3/2 , m+1/2 
\j +V] n n,j + n n,j n n~l,j + ft n-l ,j 



I m+3/2 m+3/2 
«2 I ^nj — i rn,j-l 



, m+3/2 , m+3/2 
%+l,j "I" "-raj 



(A5) 
(A6) 
(A7) 

(A8) 



• The first two equations allow for direct time advance of e n j and (p n j 

-1 



m+l 



, L + " 

l d T 2 

m+3/2 _ m+1/2 



e n ,j - 9] sm (p n<j + 



, m+1/2 , m+1/2 
m+1/2 n n,j + l ~ "nj 



+ 3 / 2 „r,A um+3/2 



m+l 



, m+3/2 



• Substitution of h n ■ and h n _ X A from Eq. (A8) into Eq. (A7) leads to the tridiagonal linear system for 



m+3/2 m+3/2 m+3/2 m+3/2 

, m+1/2 7 m+1/2 

fc n,, + * — " ~ A " — + Kj -K-IJ 



for n = 2, . . . , N with k™^ 3 ^ 2 = 0; k^^j = 0. Solving this system, we advance fc n ,j 



■^m+3/2 
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After finding fc™^ 3 ^ 2 , we update h™j 3 ^ 2 using Eq. (A8). 



The long-range instabilities are only expected for very 
tall stacks N > 1000 which are very difficult to simulate 
directly. To probe these instabilities, we use the coarse- 
grained model. Assuming that the perturbations are 
smooth in z-direction we introduce a discretization step 
Sz containing many junctions, Sz S> 1, and write coarse- 
grained equations only for junctions with n = Szm, 







dr 


= -v c e^ 








dr 




dk m 


l 


dr 


Vab 







v c e m - g(u) sin</p m + 



dh m 
du ' 



h m h m —i 
Sz 

f dfm k m+ i — k 
Sz 



Transforming variables, u = Szu, h m = h m Sz, we ar- 
rive to original equations ( Al l-( |A4| with replacements 
u u, h rn —} h m , and I — > £ cg — IjSz — X/Szs. There- 
fore, simulations of the same model with the smaller pa- 
rameter £ is equivalent to coarse-graining and allows us 
to explore the long-range instabilities in very tall stacks. 
Renormalization implies that the stack width is now mea- 
sured in units XjSz. The effective stack height is given 
by N = SzN cg , where -/V cg is the total number of the c- 
axis slices in the coarse-grained model. Having in mind 
to probe the long-range stability of the kink state, we 
use the modulation function g(u) = sgn(u — L x /2). With 
such modulation function the system is stable with re- 
spect to the short-scale perturbations. 
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